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CHAPTER I 


INTRODUCTION 

It has been recognized for some time that proper calculation of 
supersonic internal and external flows requires consideration of the 
mutual' Interaction between the highly viscous region near the surface 
and the weakly viscous region away from the surface. For supersonic 
internal flows such as in mixed compression inlets or in exit nozzles 
it is furthermore desirable to be able to track the shock structure 
even into the boundary layer to properly incorporate design features 
such as for example boundary layer control. 

Over the past decade much interest and effort have been devoted 
to viscous-inviscici interaction analyses. However, most of these 
methods have failed adequately to predict experimental data on con- 
figurations with significant curvatures or in situations such as at 
high Mach numbers where the boundary layer is relatively thick. 
Generally there are two reasons for this failure: a) The interaction 

procedure between the inviscid supersonic region and the highly 
viscous layer near the wall may have some impropriety, and b) the 
description of the viscous layer in the classical boundary layer manner 
is inadequate when curvature effects are important. 

There are several approaches for attacking the viscous-inviscid 
interaction problem in supersonic flow. The ultimate approach would 
be to solve the time dependent Navier-Stokes equations for the whole 
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flow field thus eliminating the patching procedure between the viscous 
region near the surface and the inviscid region away from the wall. 
However since the Navier-Stokes equations are spatially elliptic it 
is unfeasible and impractical from the points of view of both com- 
puter storage and running time. The objective of the present work 
is to develop a viscous-inviscid interactive procedure in supersonic 
flow that represents an intermediate development between past treat- 
ments and exact solution of the Navier-Stokes equations. The inter- 
active flow analysis developed herein is based on dividing the flow 
field into regions where in one region the flow is supersonic with 
a dominant inviscid character and is treated by a hyperbolic system 
of equations , while the second region where the flow is highly vis- 
cous is treated by a parabolic set of equations. In both regions 
forward marching techniques can be used thus considerably reducing 
storage and time requirements. However when replacing the Navier- 
Stokes equations that are elliptic with sets of equations that are 
hyperbolic and parabolic in character, there is no capability of 
directly dealing with upstream influence effects. Hence the present 
procedure cannot handle flow separation or strong shock interactions 
where the elliptic behavior is important. 

A. Applicable Earlier Work 

The interaction methods that have been developed for supersonic 
flow all consider the flow field to consist of two regions. Generally 
these regions are a boundary layer or equivalent viscous layer and 
an external supersonic region. The methods however differ in the 
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complexity of the modeling assumptions for each of the regions and 
in the coupling procedures for effecting interaction. 

An early but most significant formulation of an interaction 
procedure for supersonic flows was by Crocco and Lees (reference 1) . 

In their work they related the pressure distribution of the external 
supersonic flow to the local slope of the displacement thickness of 
the viscous region using the Prandtl-Meyer relation which is a one- 
family characteristics procedure. The viscous layer in the Crocco- 
Lees procedure is based on the classical boundary layer approximation 
but was treated as a mixing layer using a specially-developed momentum 
integral procedure. Lees and Reeves (reference 2) extended the Crocco- 
Lees method by additionally employing a moment of momentum integral 
equation to improve the treatment of entrainment. The extension of 
the Lees-Reeves integral interaction procedure to include consideration 
of heat transfer was by Klineberg and Lees (reference 3) . 

Reyhner and Flugge-Lotz (reference 4) improved the treatment of 
the viscous portion of the interaction analysis by applying a full 
finite difference technique to solution of the compressible laminar 
boundary layer equations in the physical plane. The boundary layer 
is treated in the classical limit with the normal pressure gradient 
taken as zero. As with the earlier described procedures, the coupling 
between pressure and local streamline deflection at the edge of the 
boundary layer is through the Prandtl-Meyer relation. 

Miller (reference 5) argues however that the inviscid flow must 
be calculated by a two-family characteristics method in order to 
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obtain a mathematically well-posed supersonic interactive problem, 
and thus eliminate the saddle point type singularity that is intro- 
duced by using the Prandtl-Meyer relation which is a one-family 
characteristics solution. 

Ferri and Dash (reference 6) improved the treatment of viscous- 
inviscid interactions in supersonic flow in two ways: first, by 

applying a higher approximation for the boundary layer that includes 
normal pressure gradient and longitudinal curvature effects, and 
second, by using a two-family rotational characteristics scheme in 
the outer region that allows for entropy changes due to viscous effects. 
The pressure distribution across the viscous region was assumed to 
be a polynomial of fourth degree uncoupled from the rest of the system. 
The coefficients were determined by assuming that the first and 
second normal derivatives of the pressure at the wall are zero and 
that the remaining terms are de,>endent on the longitudinal curvature 
effects. The system of equations obtained were solved numerically. 

In the viscous region the x- momentum and the energy equations were 
expressed in finite difference form and solved simultaneously for u 
and T. The normal velocity distribution was obtained by integration 
of the continuity equation and the process repeated iteratively until 
convergence was obtained for u, v, T and p. 

B. The Present Method 

The present work is an extension of the idea of Ferri and Dash 
(reference 6) wherein the flow field is divided into two regions: 
a) an inner region which is highly viscous and mostly subsonic and 
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b) an outer region where the flow is supersonic and in which viscous 
effects are small but not negligible. 

The inner region is treated by a system of equations of the 
boundary layer type. This system is obtained by reexamining the 
Navier- Stokes equations for steady compressible two-dimensional or 
axisymmetric flow in curvilinear coordinates and through an ordering 
procedure retaining terms of order unity and (6/L). In addition to 
the classical boundary layer equations the system of equations so 
obtained includes in a consistent way the second order effects of 
longitudinal and transverse curvature as well as normal pressure 
gradient . 

In this system the normal momentum equation is retained. The 
equations are a coupled parabolic set in the longitudinal velocity, u, 
the normal velocity component, v, and the static temperature, T. By 
incorporating a suitable effective viscosity hypothesis, the system 
can be used to calculate both laminar and turbulent boundary layers. 

The system of equations obtained is solved simultaneously in the 
physical coordinate plane using an implicit finite difference technique. 
This procedure provides an exact and stable numerical solution to the 
viscous flow equations in the inner region. 

The numerical solutions for the outer region are obtained by 
applying the method of characteristics to a system of equations 
which includes viscous and conductive transport terms normal to 
streamlines. In this streamline-normal coordinate system, terms 
of order unity and (6/L) are retained for the viscous and heat flux 
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terms added, whereas curvature effects are kept fully. By introducifig 
the transport terms as corrections, the equations retain their hyper- 
bolic character. These correction terms include additional second 
order terms, over and above those retained by Ferri and Dash (reference 
6) . The solution of the characteristic equations have been structured 
as an inverse grid scheme in a streamline-normal network. In this 
reverse scheme both characteristic Mach lines are extended back until 
they intersect the known data region on a normal to the streamlines. 

The streamline condition has been replaced by a streamfunction con- 
dition thereby preserving mass flow within a stream tube. This allows 
for a very equitable mesh distribution which always maintains itself 
in the downstream direction without redistribution of the grid points. 
The resulting system of equations in both the outer and the inner 
regions are consistent to order (<S/L). 

In the interactive mode following the suggestion of Ferri and 
Dash (reference 6) , the inner and the outer regions are matched 
along a line where the Mach number is approximately 1.2. The match- 
ing conditions are continuity of the flow variables u, v, T and p at 
the interface. The detailed algorithm of the interactive procedure 
for the interaction mode is given. 

Each of the portions of this analysis will be discussed separately . 
The development of equations and numerical solution procedure for the 
outer region together with some illustrative examples is presented 

j in Chapter II. This is followed in Chapter III by an equivalent 

1 presentation for the inner region. The interaction procedure between 
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the two regions is described in Chapter IV, Discussion and summary 
of the major portion of the present work is given in Chapter V. 


CHAPTER II 


OUTER REGION 


A. Equations of Motion 

The equations of motion for the outer region are written for 
steady viscous compressible two-dimensional or axisymmetric flow. 
These equations for the coordinate system shown in figure 1 are: 


Continuity 


3 (p*u*y*'^) ,3 (P*v*y*'^) 

+ - 0 


3x* 3y* 

Longitudinal momentum 


^ . 9u* , . . 3u* 

p*u* + P*v* 


3x* 


3y^ 


1 


3x* y*^ 3y* 


3 / jP 9u*. 


3y*' 


+ A (p* -1^) 


3x* 


3x*' 


Normal- momentum 


p*u* 


9v* . j. , 3v* 
-3^ 


+ p*v* 
3 

3x* 


3p* 

3y* ' y*^ 9y* 


+ 


3 / ju PI.CF 3v*,, 

(h*y* 


3y' 


(p* — ) 
3x*^ 


Energy 


p*u*c* + p*v*C* - (u* + V* 

^ p 3x* p 3y* 3x* 3y* 

1 3 , .0, . 9T*, . 3 . 3T*, 

(y* k* — ) + (k* ~) 


y*n 9y* 


3y*' 


3x* 


3x*^ 


^ ^3y* 3x*^ ^ ‘•'^3x*^ ''3y*'^ ^ 


(1) 


( 2 ) 


(3) 


( 4 ) 
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Egtiation of State 

p* = p*R*T* (5) 

In these equations 

0=1 for axisymmetric flow 
0=0 for two-dimensional flow. 

The equations of motions (1) to (5) are made dimensionless in 
the following manner: 


u^ 


y = 


\et 


V = 


V* 


u 


■; p = 


_ 21 


REF 


REF 


P = 


T = 


X* 


; y 


c =-^ 
^REF P ^ 


k = 


REF 

k* 


REF 


X = 


X* 


%EF 


REF 


REF 


( 6 ) 


The resulting dimensionless equations are: 
Continuity 

3 (puy*^) a_ (pvy^) 


3x 9y 

Longitudinal-momentum 
9u , 3u 

p u + p V— 

3x y 
Normal- momentum 


= 0 


(7) 
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■1 3 


I^EF 


'REF 
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3x^ 3x^ 
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3v , 3v 

3x " 


lE 


• 1 3 




o 3v > 
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{-- — (yy^ — ) + } 

^ 3y^^^ 3y^ 3x^^^ 
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Energy, 


puC 


8T 
p 3x 


_ 3T (y- 1) , 3p , ^P^ 

+ pvC r ^(u ^ + V -^) 

p 9y Y dx 9y 


"Pt* H.0 

REF REF y 


o 3y ' 3y Pr 


Re 

REF REF 


9x 9x 


REF 


( 10 ) 


Equation of state 
p = pT 

where: 


\ef 


REF 




1/2 


Re 


REF 


Fr. 


= Prer^ref\ef 

*^REF 


C 

Fref 


Reference Mach number 


Reference Reynolds number 


Reference Prandtl number 


(11) 


'S(EF 

The equations of motions (7) to (11) are now transformed from cartesian 

coordinates (x,y) to curvilinear coordinates (s,n) where: s and n are 

respectively the distance along a streamline and the distance normal 

to a streamline. Transport effects such as viscous shear and heat 

flux are retained only normal to the streamline. This is because from 

92^2 2 

an ordering procedure, terms such as Tr~vl~7r~2. order of (6/L) 

a s on 


that can be neglected in the present analysis. 
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The governing equations (7) to (11) when transformed and simplified 
(details given in Appendix A) are: 

Continuity 


9 / \ , 90 _pq . Q 

■— (pq) + pq T- = - o „ sxne 
9s 9n y 


s-momentum 


9q + — ^ - Q, 
pq -.,2 9s ] 

^“ref 


n- momentum 

2 96 ^ 

ir-*- 


Energy equation 


.2 9n ^2 


pq 


H _ ix-ll , lE = 


9s 




(12) 


(13) 


(14) 


(15) 


Equation of state 


p = pT (16) 

where q is the velocity in the streamline direction, 6 is the stream- 
line direction and Q^, Q^, are correction terms due to viscous 
shear and heat flux. The detailed expressions for the correction 
terms are rewritten here: 


1 ^ MOS0 9£ - pq(-|^) - sin^e |^} 




(17) 


‘2 Re 


REF 


r 9q 99 , ^ dO , d , dD.-. 

in * y ® ’ 


3. li + ii; 

9n 9n^^^ 3n' 


(18) 




( 19 ) 
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By considering these terms as known source terms the systems of 
equations (12) to (16) retains its hyperbolic character. 


B. Method of Characteristics 

The characteristics derived from the system of equations (12) 
to (16) are defined by two ordinary first order differential equations 
where the independent variables are the static pressure, p and the 
streamline direction, 6. The detailed derivation of the characteristics 
relations is given in Appendix B. 

The characteristics equation is: 

= ±tan\ (20) 


where X = sin ” is the Mach angle. These characteristic directions 
M 

are the same as for inviscid flow. The compatability relation is 
given by: 


iP- 

yP 


d0 ^ , osin9 *^1 

sinXcosX y 2 

pq 


(y-i)m 

pq 


REF 




Q_M sinX 
2 REF 


dx 

cosXcos (6±X) 

( 21 ) 


The + and - signs correspond to C and C characteristic lines res- 
pectively. The quantities Q^, Q^j functions of the local normal 

gradients of velocity and temperature and have the same role as the 
entropy terms in equivalent equations for inviscid rotational flows. 

The variations of entropy, S , and stagnation temperature along 
the streamline as derived in Appendix B are: 


= 0 


d S ^ ~^\eF 
ds pq 3 


( 22 ) 
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ds pq ^3 p 


(23) 


The oblique shock wave relation required in dealing with shock 
points in the characteristic net is obtained from Reference (7) and is 
given by: 


6 4 2 
sin ^ ^SH ^ ^SH 


= 0 


(24) 


where 


3 is the shock angle 
SH 

, M^+2 . 2 . 

b = - — ^ ysxn o 

2M^+1 , t/Y+In^. (Y- l)n . 2. 
c = ^ + [ (-2 ^ ^Jsin 5 

M M 


d = - 


cos 6 


M 


In these equations, M is the Mach number upstream of shock wave and 
6 is the local streamline deflection angle. 

In order to establish a well-posed problem, the following boundary 
conditions must be given: 

(a) data for all quantities must be prescribed along 

an initial datum line. In the absence of the 
correction terms Q-i »Q2 Q^, only p, 6 and M need, 

be specified to allow the characteristics calculation 
to proceed. However, upon including Q^, Q2 and Q 
the velocity and temperature information is needed, as 
well. 

(b) the shape of the boundary surface or bounding stream- 
line, y = Yg(x). 

C. Numerical Procedure 


The conditions at each grid point (x,y) in the physical plane 
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are determined by the characteristic equation (20) and the compatibility 
condition (21) . A non-uniform grid point distribution on the normal to 
the streamline is chosen to allow a fine mesh spacing in the inner 
portion of the outer region and coarse spacing at the outer portion 
of the outer region. The computational procedure divides itself into 
the following four basic elements: 

1. Conical flow field calculation (for axi symmetric flow) 

2. Boundary point calculation 

3. Field point calculation 

4. Shock point calculation 

Conical Flow Field Calculation - The flow past an infinite cone 

with attached shock front is described by the Taylor-Maccoll equations 

which are solved numerically by means of Runge-Kutta integration (see 

Appendix C) . The conical shock angle, is calculated by an iterative 

SH 

procedure beginning with a first guess of the shock angle for the given 
freestream conditions. The integration is continued to the cone sur- 
face and the resultant cone angle compared with the specified cone angle. 

A new estimate of the shock angle 3 is made based on the error and 

dH 

the process is repeated until convergence is obtained. The converged 
solution gives complete information on the conical flow field. 

Boundary Point Calculation - Calculation of a boundary point 
(Fig. 2) requires that the normal velocity at the boundary point be 
zero. The boundary point A can be either on a solid boundary or a 
point on a prescribed streamline whose deflection is known. The 
pressure at the boundary point A is calculated in an iterative way. 

Using a reverse scheme the C- characteristic is extended from point A 
in the upstream direction until it intersects the normal to the 
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streamlines at point B. The flow variables at point B are determined 
by interpolation. Using average values of points A and B the com- 
patibility condition for the C- characteristic curve will give new 
flow variables for point A in terms of the known values at point B. 
Iteration is continued until convergence is obtained for point A. 

Field Point Calculation - Calculation of point D (Fig. 3) is 
done using the reverse scheme, that is both the C+ characteristic 
and the C- characteristic lines are extended from point D in the up- 
stream direction where they intersect the upstream datum line at points 
A and B respectively. Point D is located on the normal to the stream- 
line through point E [x(2,J-l), y(2,J-l)]. 

The compatibility equations at point D determine the flow 
variables and the streamline slope for point D in terms of the flow 
variables at points A and B. Integrating the mass flow between point 
D and the previously calculated point E [x(2,J-l), y(2,J-l)] the 
stream function at point D is: 

In order to preserve mass, the new location of point D on the normal 
to the streamline, is corrected with respect to the difference de- 
fined by: 

6ijj = - ip(J) (26) 

Normally convergence is obtained in a few iterations. 

Shock-Point Calculation - Shock point calculations are basically 
different than the field point calculations. Since conditions upstream 
of the shock wave are known, the oblique shock relations must be in- 
corporated into the calculation procedure. 


Point B (Fig. 4) is located just downstream of the shock wave 
and on the normal to the streamline from point C[x(2,J-l), y(2,J-l)]. 

An initial guess is required for the flow variables and streamline 
deflection at point B. Then point A can be determined by the inter- 
section of the C+ characteristic through point B and the upstream normal 
to the streamlines. The shock wave angle is calculated using equation 
(24). This determines the flow variables at point B from oblique 
shock wave relations. The coefficients of the compatibility equation 
are updated using averages of the flow variables at points A and B. 

The streamline deflection is recalculated and thus a new shock wave 
angle is calculated. Again in an iterative way, convergence is 
obtained. 

D. Numerical Results 

To test and illustrate the procedure developed for the outer 
region, two sets of calculations were performed. The first of these 
is for the waisted body tested by Winter, Smith and Rotta (reference 
8) while the second set is for a Mach number 3.5 mixed compression 
inlet for which characteristics calculations were performed by 
Syberg and Hickcox (reference 9) using a different scheme than 
developed herein. 

The geometry of the waisted body of revolution is described 
(figure 5) by a set of five polynomial functions each pertaining to 
only a section of the body. Calculations were performed for a 
series of supersonic Mach numbers for which experimental data are 
available. The correction terms Q^, and have been excluded 


17 


in these calculations. 

Figure 6 shows the oblique shock waves for M =1.7, 2.0, 2.4 

00 

and 2.8. The wall Mach number distributions along the body are shown 

in figure 7 for the same free stream Mach numbers. A comparison of 

the calculated Mach number distribution on the surface and the measured 

Mach number distribution at the edge of the boundary layer (reference 8) 

shown in figure 7 shows very good agreement except perhaps at = 2.8 

in the vicinity of — a? 0.4. Examination of the experimental 

^REF 

velocity profile at this Mach number and location (figure 9 of reference 

8) shows that the velocity is still slowly increasing beyond the nominal 

"edge” of the boundary layer as chosen by Winter et al (reference 8) . 

The calculated wall static pressure distributions for the same 

free stream Mach numbers are shown in figure 8. It is seen that the 

pressure minimum occurs consistently in the vicinity of the inflection 

point of the body (- :^i0.42). Unfortunately, experimental wall 

^REF 

pressure data are not available in reference 8 for comparison with 
these results. 

The streamline patterns as calculated for the waisted body at 
M^ = 1.7 and M^ = 2.8 are shown in figures 9 and 10 respectively. Of 
parti c:utb.r interest is the appearance at = 2.8 of a noticeable 
second shock wave in the flow field emanating from x/L 0.75. 

ivEr 

This wave for the same Mach 2.8 test of Winter, Smith and Rotta was 
also identified in the "smooth shock-capturing" calculation of 


reference 10. 
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Figure 11 shows the pressure distribution along normals to the 
streamlines emanating from a number of stations on the body. The 
presentation for = 1.7 (figure 11a) is a composite of the pressures 
along the normals to the streamlines labeled A through N in figure 9. 
The ordinate in figure 11a is (y-y=,-J / (y„,.-y„) for the points on the 

15 bn. B 

aforementioned normal lines. From the pressure distributions for 

X /L _ > 0.787, a mild shock wave seems to appear in the pressure 
B REF ~ 

distribution although there is no apparent corresponding streamline 
deflection in figure 9. 

Figure lib gives a similar portrayal for = 2,8. The shock 
wave appearing in the pressure distribution for t 0.787 is 

readily identified with the locus of streamline deflections in figure 
10 and is coincident with the shock location as calculated in ref- 
erence 10. 

Numerical calculations were also done for the Mach number 3.5 
mixed compression inlet sketched in figure 12. Geometrical data for 
this inlet are given in reference 9. The pressure distribution 
calculated with the inlet by both the present streamline-normal 
procedure excluding the correction terms Q^, and and by the 
more conventional characteristics scheme of reference 9 are compared 
in figure 13 for the center body surface and in figure 14 for the 
cowl surface. The excellent agreement between the two computational 
schemes indicates that little if any accuracy is lost by neglecting 
Mach wave intersections between the datum planes that is characteristic 
of inverse scheme methods. 
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In order to investigate the effects of the correction terms in 
the outer region, a power law velocity profile (n=ll) of suitable 
thickness has been imposed in the Initial datum line where the wall 
Mach number was chosen as Mw = 1.9. The wave pattern for the mixed 
compression inlet is calculated both with and without the correction 
terms. For the purposes of this exercise, the correction terms 
are evaluated using Sutherland viscosity and Pr = 0.72. The 
numerical results shown in figure 15 indicate a shift of the wave 
location and curving of the wave reflection due to the viscous and 
conduction terms. This partial simulation of the presence of a 
boundary layer gives some indication of what might be expected in an 


interactive calculation. 


CHAPTER III 


INNER REGION 

The various methods for dealing with the viscous region differ 

in the order of the terms kept in the system of governing equations. 

In past work on interactions the tendency has been to use classical 

boundary layer equations (with ■” = 0) or else such equations augmented 

dy 

by a centrifugal correction for normal pressure gradients. The pre- 
sent interaction procedure depends however on having accurate re- 
presentations of the normal velocity distributions at the matching 
location in the boundary layer even in situations with longitudinal 
and transverse curvature. Hence it is desirable herein to employ a 
set of equations that are consistent to order (6/L) relative tb the 
leading or classical boundary layer equations. In that way the sets 
of equations for the inner and outer regions are consistent to order 
(6/L). Within this second order approximation, the normal momentum 
equation appears as a coupled member of the set of inner layer equations 
Before proceeding into the details, a brief review of pertinent second- 
order effects is appropriate. 

The bo undary’’ layer concept introduced by Prandtl has been 
successful in yielding solutions for viscous flows at high Reynolds 
numbers, as long as the boundary layer remains attached to the surface 
and remains thin enough so that it does not noticably affect the 
external flow. A general discussion of higher-order approximation 
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for boundary layers has been given by Van Dyke (reference • 11) . In 
that review article, Van Dyke classified the second-order corrections 
into two categories according to whether the additional terms appear 
in the differential equation through the curvatures or through inter- 
action with the external flow. Van Dyke further subdivided each 
category giving some analytical details and physical interpertation 
for the described effects. 

Maslen (reference 12) has presented a set of equations that is 

complete with respect to the inclusion of second order boundary- 

layer effects. This set of equations is the starting set for the 

developments in this chapter. Maslen himself in reference 12 went 

on to study some weak interaction questions by means of similarity 

solutions. Seginer (reference 13) solved a system of equations 

. 1/2 


retaining terms to order (• 


-) where Re 


^REF 


is the reference 


length Reynolds number. In his system the normal momentum equation 
is a coupled member of the parabolic set of governing equations. 
Seginer then went on using the Stewartson transformation and 
similarity arguments to obtain ordinary differential equations for 
the second order boundary layer system. He obtained solutions for 
a two-dimensional flat plate at Mach number 4 to Illustrate the 
very slight effect even for tbat case of the normal pressure gradient. 

The present system draws on both of the above works in its 


development of the second order system of equations but the solution 
procedure is developed in physical coordinates for later interactive 
matching with the outer solution. Similarity arguments are notitivoked. 
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A. Equations of Motion for the Inner Region - Laminar Flow 

The highly viscous flow in the inner region is assumed to be 
steady, compressible, laminar or turbulent, two dimensional or axisymmetric 
flow over an adiabatic or non- adiabatic surface. The appropriate 
system of equations in dimensional form is given in reference 11. 

Expressed in curvilinear coordinates (?,?) in which 5 is measured 
along the surface and ^ normal to the surface (figure 16) and where 
u, V are the corresponding velocities, these equations are: 

Continuity 

+3^Jp*v*(l +g-)l - 0 (27) 
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g^-momentum 
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( 30 ) 


In these equations r* is' the lateral radius of curvature of the point 
(5*,g*), R* is the longitudinal radius of curvature of the surface, and 
is positive for convex surfaces (reference 11) while 
cf = 0 for two dimensional flow 

and 

a = 1 for axisyminetric flow. 


4 - 
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The angle 6 is the local slope of the axisyiranetric body relative to 
w 

the normal to the axis. Thus 

cose = 

w d?* 

Reexamining the equations of motion using dimensional analysis 
(see Appendix D) , a set of equations is obtained constituting a higher 
approximation than the classical boundary layer for the inner region, 
since terms of order unity and of order (6/L^^p) are retained. The 
non-dimensional equations of motion obtained to this order expressed 
in terms of reference Reynolds number become: 

continuity 


|j[pur“] +7T72 R>) = “ 


Re 


REF 


(31) 


g- moment urn 
3u 


M 4. - ^ 4 . ^ -BdiX 4 . i£ 

pu — + pv(l + jj) + E + 

^ Re 


Re 


REF 




REF 


13 , O'/I _L 

a 34^^’' „ 1/2 R' 3? 

r Re 


1 is Ju, 

/ 'S J 


1 u 3 , Ov 


'REF 


„ 1/2 a R 3? 


(32) 


t-momentum 


1 r 9v , . 1 pu j_ t,„1/2m X 1 

9g 1/2 R^ 3?^ ■ R ^ REF^^ 1/2 R^ 3? 

^^REF ^^REF ^^REF 


2 13/ o iy, . 1 _L 3_, , o ^ 

1/2 a 3?^^’' 3C^ „ 1/2 a 3g^^'' 3?^ 

Re r Re„„„ r 


'REF 


'REF 


2 l/2^r 3u^ 3v , , 

3 ^®REF 3C*- ^ 3g 3C r ^ 


( 33 ) 


25 


energy 
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B . Equations of Motion for the Inner Region - Turbulent Flow 

When the equations of motion are written in terras of mean quantities 
(velocities, pressure, density, temperature, etc.) and fluctuations 
about the mean, and then averaged with respect to time, the resulting 
set for turbulent flow (a detailed development is given in Appendix E) 
keeping only the leading correlations u'v* and v’T' is: 
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C. Viscosity Laws 

The equations are formulated to accommodate any variation of 
viscosity that may he required in properly implementing a boundary 
layer calculation. For the examples presented in the work the dynamic 
or absolute viscosity is represented by the Sutherland viscosity 
relation which in dimensionless form is 
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y = 




(41) 


where: 

198.6 

a = 

T 

REF 

For turbulent boundary layers, the momentum and energy transfer 
correlations u’v' and v'T' are included in the equations of motion and 
are evaluated using Boussinesq scalar eddy viscosity and eddy con- 
ductivity coefficients through which these transport terms are in 
turn related to local mean velocity and temperature gradients. Thus 
the turbulent transport coefficients are defined: 


-v’T’ = 


3u 
e . ~ 
u 3C 

^t 3C Pr^ 9C 


(42) 


where e is the scalar eddy viscosity and Pr is the turbulent Prandtl 
u t 

number . 

Eddy Viscosity Model : Experimental data with equilibrium tur- 

bulent boundary layers indicates that the scalar eddy viscosity 
function can be simulated by a two-layer model (reference 14) . The 
inner layer in the vicinity of the wall is characterized by increasing 
turbulence with distance from the wall namely varies almost 
linearly* with distance from the wall. In the outer layer the scalar 
eddy viscosity function is nearly constant except for the intermittency 

factor^ Cebeci and Smith (reference 15) have extended the Van-driest 

______ ^ 

In the very near neighborhood of the wall 'u y . 
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formulation of the law-of-the~wall (reference 16) in order to include 
effects of pressure gradient and heat and mass transfer at the wall. 

The Cebeci-Smith model for effective eddy viscosity will be used in 
the present calculations for turbulent flow. While these are other 
suggested models for the eddy viscosity (reference 17 for example) , 
the Cebeci-Smith formulation has been shown to be adequate for engineer- 
ing calculations. 

The Cebeci-Smith model (reference 15) written for axisymmetric 
flow* is as follows : 
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where : 
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A-is a function of pressure gradient, mass transfer at the wall and 


viscous shear stress at the wall and is given by: 
w w 

where 
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The two-dimensional version can be extracted by settxng — = 1 + 
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and evaluating the expressions in the limit as r 
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Rg = Reynolds number based on momentum thickness 
The intermittency factor is given by: 
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The condition of continuity of the eddy viscosity function determines 


the values of , i.e. : 
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( 46 ) 


The two- layer representation of the eddy viscosity function is 


sketched in figure 17. 


Substituting the expressions (42) for the turbulent transport 
terms into equations (37) to (39) the system of equations becomes: 
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equation of state 


P = 






(51) 


Equations (31) to (35) for the laminar case or equations (47) to 
(51) for the turbulent case are a Cvupled parabolic system in the 
variables u, v, and T, hence two pointy boundary conditions for u, v, 
and T and one boundary condition for the pressure are required. For 
convenience in computation the pressure is split as follows: 


p(e,C) = + Pt(?.?) 


(52) 


EXT 

where p„^^(C) is the external pressure (presumably known), and 
EXT 

Pj(?>C) is an induced pressure due to normal momentum consideration. 


D. Boundary Conditions 

The boundary conditions at the wall (? = 0) are: 


u(^,0) = 0 
v(C,0) = 


for impermeable walls 
^w^^^ for suction or blowing 


T(5,0) = T (C) 
w 


or 


ill 


= 0 


(53) 


The outer boundary conditions depend on the use made of the inner- 
region equations. Those presented here are the outer conditions for 
boundary layer calculations. The outer boundary conditions employed 
in the interaction procedure are described in Chapter IV. For com- 
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putational purposes, the outer edge of the boundary layer is taken to 
be well outside the conventionally defined thicknesses. In practice 
this turns out to be three to four times 6* for the laminar boundary 
layer and of the order 106* for the turbulent boundary layers. At 
this location: (S) 


Pj(C,<S) = 0 


Reil: ^ C-6 


^ 6^5 d? 
dT, 


d? 


d? ^6 d^ 


REF 
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Initial Profile : For starting the finite difference calculation, 

profiles for the unknowns, u, v, T and p at a specified initial 
station are needed. For laminar boundary layers, similarity 
solutions have been incorporated to start the finite difference flow 
field calculation. For incompressible turbulent boundary layers the 
starting profile has been constructed from the law-of-the-wall. 


E. Numerical Solution of the Inner-Region Equations 

For solution of the non-linear partial differential system, 
equations (31) to (35) for laminar flow or equations (47) to (51) 
for turbulent flow are linearized and then replaced by a system of 
linear algebraic equations using a modified Crank-Nicolson implicit 
finite difference scheme. Since because of the nonlinearity of the basic 
system the flow variables in the coefficient matrix depend on the 
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solution vector, an iterative procedure is applied until the differences 
between the flow variables for successive iterations is as small as 
desired. 

A primary objective in the development of a numerical procedure 
is to get it to yield a stable and convergent solution for the system 
of finite difference equations. Stability and convergence of 
numerical solutions of partial differential equations is discussed by 
Roache (reference 18). Basically, instability results from un- 
avoidable small perturbations in the flow field due for example to 
round-off error, truncation error, etc. If in marching downstream 
the errors diminish then the method is stable; if the errors grow in 
inarching downstream the method is unstable. 

It was decided to apply an implicit difference schema rather than 
an explicit one. A broad discussion of the two schemes is given in 
reference 18. Generally it is expected that implicit difference 
schemes are far more stable than explicit difference schemes. On 
the other hand implicit schemes involve a system of algebraic 
equations that must be solved simultaneously since the equations 
are coupled. 

Other investigators have attempted solving the flow equatxons in 
the physical plane using implicit difference schemes. Reyhner and 
Flugge-Lotz (reference 4) solved the system of classical boundary 
layer equations simultaneously in the physical plane for u, v and T, 
Their numerical results indicated an oscillation in the normal 
velocity at the edge of the boundary layer. Such oscillations are 
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undesirable if the method is to be used as part of an interaction 
calculation. Since their system of equations includes only continuity, 
x-momentum and energy equations the normal pressure gradient is 
assumed to be zero. 

Ferri and Dash (reference 6) have included an uncoupled normal 
momentum equation in their mathematical model. Their numerical 
solution was also done in the physical plane. Although they used an 
implicit difference scheme, their procedure is quite different than 
that of Reyhner and Flugge-Lotz (reference 4) . They assumed the 
pressure distribution across the boundary layer to be a polynomial 
of fourth order uncoupled from the rest of the system. The co- 
efficients were determined by assuming that first and second normal 
derivatives of the pressure at the wall are zero and that the re- 
maining terms are dependent on the longitudinal curvature effects. 

The system of equations i.e., continuity, x-momentum and the energy 
equations were expressed in finite difference form. Solution of the 
finite difference equations was done successively; that is they 
first solved simultaneously for u and T using the x-momentum and 
energy equations, then the normal velocity distribution was ob- 
tained by integration of the continuity equation. The process 
was repeated iteratively until convergence was obtained for u, v 
and T. The converged solution was then used to determine the variation 
of p across the boundary layer. 
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The objective in the present work is to solve the second order 
system of equations for the viscous region in the physical plane and 
by extending the Reyhner-Flugge-Lotz method (reference 4) and 
solving the coupled parabolic system of equations simultaneously. 


Difference Scheme and Quotients - The grid scheme on which the 
expressions for the difference quotients are based is shown in 
figure 18. A grid with variable mesh size in the normal 4 direction 
and a uniform mesh size in the longitudinal ^ direction has been 
chosen. In order to obtain a fine mesh near the surface where the 
gradients are large and a coarse mesh away from the surface, a 
geometric series has been chosen to locate the grid points in the 
c. direction. 
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K - 1.0 
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(55) 


where: 

DK = first interval 

K = ratio between two consecutive intervals 
= C coordinate of the n^^ grid point 

The truncation errors of the difference quotients are based on 
a Taylor series expansion of a function with two Independent vari- 
ables about a point where the function and its derivatives are known. 
The validity of the expansion depends on the existence and continuity 
of all derivatives of the function f(x,y) which here represents any 
of the unknown functions u, v, T and p. The value of f(x+-h,y+5,) 
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can be expressed as follows: 
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For any of the unknown functions u, v, T and p represented here either 

by f(x,y) or g(x,y) the expressions for first and second derivatives 

of the above mentioned functions at the point A (figure 18) located 

at (? , ) depend on the differencing scheme. The expressions 

nrfl/2 n 

for centered differences, backward differences and forward differences 
are given below: 

Centered Differences 
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subtracting (58) from (57) yields 
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2) First derivative in ^-direction: 

The derivative in ^-direction is obtained as a weighted average 
of the derivatives in the ^-direction at station m and station m+1. 

+1 = + If I + 7T Ml " 

mfl,n+l nrfl.n 3^' 2! n+1 

’ ’ m+l.n 3? m+l,n 


= f _ If I Ar + -i- ^^^ 1 

nrl-l,n-l m+l,n 3C . . n 2! ^„2' 

m+l,n 3C nrt-l,n 


(AC f + 
n 


Since Ac = Ac K, the resulting expression is 
n+1 n 


nrfl,n “'=’n 

Similarly for station m 


f f 

m+l,n+l m+l,n-l 

AC^a+K) 


(?[(AC ^)] 


^m,n+l ^m,n-ll 


AC (1+K) 
n 


Finally we obtain: 


Ifl = K (ff) I + (1-X ) (|f) I (64: 

^ m4-l/2,n ^ nri-l,n nijii 

For = 1/2, the Crank-Nicolson centered difference expression is 

obtained. 


Second Derivative - Second derivatives appear in the flow 
equations only in the C-direction. For centered differences: 


^nri-l,n+l ^m+l,n ^^n+1 2! 

(65) 


. ^ 9f| * . 8 f| n2 

^nrfl,n-l ^m+l,n ~ 3cl ^n 2! .^ 2 ! ^ 

’ ’ nrfl,n 3C nrfl,n 

( 66 ) 
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Multiplying equation (65) by (AC^) and equation (66) by (AC^^j^) and 


adding the resulting equations yields 


3? mfl,n 


nrfl.n+l ttrH,n-l m+l,n 

K(1+K)(AS )^ 
n 


For flexibility a weighted average has been used to obtain the second 


derivative at point A: 


3? nrM/2,n 


2 2 
3 fi . n T 3 

'2 ^^" 2^ ~~2> 

3? m+l,n 3? m,n 


Again for the equal-weighted Crank-Nicolson result is obtained. 


Other Derivatives 


(^)l 

^3C 'nrbl/2,n 


^^iD-H,n+l ^mfl,n+l^ ^^m,n+l ^m,n-l^ 

(l+K)^(Ai;^)^ 


(f -f We -e ) 

Jll£| = nrfl.n+1 nrfl .n-1^ ^^m.n+1 ^m.n-1^ 

mfl/2,n 2(A?^)^(1+K)^ 


^^nH-l.n+l ^nrfl.n-l^ ^^m.n+l ^m,n-l^ 


2 (AC )^(1+K)^ 
n 


To be noted is the linearity of these expressions in the quantities 
evaluated at station mhl. 


Backward Differences 


First Derivative 


W,n-1 ' ' if' "'n ^ TT 7^1 „ 

’ m+l,n 3C irl-l,n 
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nrfl.n-2 m+l,n n-1 n 2! n-1 

(72) 

2 2 
Multiply equation (71) by (Ai; + A^ ) and equation (72) by (AC ) . 

n n— ± n 

Then subtraction of (72) from (71) yields: 


li 

3C 


I " n g +ZK ). . _ J IMl 

m+l,n ’ (1+K)AC^ nrfl,n-l AC^ m4-l,n-2 


+ f 


K 


n 


0 "[(A? )^1 

n 


(73) 


Hence the first derivative for point A (figure 18) using backward 
differencing and weighted averaging becomes: 


fl = S If I 

^nrfl/2,n m4-l,n 


+ ( 1 -X 3 ) Ifl 


(74) 


m,n 


Forward Differences 


First Derivative 

. 9f| , .1 

■p = P -p — At -p — — • ' * 

m4-l,n+l nr+-l,n i n+1 2! ^,,2 , 

’ ’ ^ m+l,n 9C nrfl,n 




£ 1 3f I . _L A \ I 1 ,3 f I 

m+l,n+2 m4-l,n ^n+2 ^n+1 2! 

’ ’ ^m+l,n 9C m+l,n 


(75) 

^^n+l'^^^n+2^ ^ 


(76) 


Multiplying equation (75) by (A? _ H- A? ) and equation (76) by 

n+1 n+2 

2 

(A? ,-) then subtracting equation (76) from equation (75) and ex- 
n+1 

pressing in terras of A? we obtain: 

n 


il 

9C 
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I = - ( 2 +K) ( 1 -t-K) 1 

nH-l,n (1+K)KAC^ ^2 :n4-l,n+l " nH-l,nf2 ^2 

n ' ^ n 

(77) 


Then as for backward differencing, the expression for the first de- 
rivative at point A is: 


•^1 ~ X + n-x ) ^ 

m+l/2,n nrfl,n 


(78) 


m,n 


Difference Equations 

Following the suggestion of Reyhner and Flugge-Lotz (reference 4) 
the system of non linear partial differential equations has been 
linearized in the following way: 
continuity 


+ f^[r“pv(l + -4 t5 

^^REF 


momentum 


(i) 


(pu)«|f+ 

35 35 Ke^EP 

2 

= _ lE a. r^i j. 4- £^M(i) 


3? 




Re 


REF 


3? 


+ f— — Cr^ (1 4- — - — -^) u (1 + -^) ) 1 


^®REF 


(i) 


"=REF 


(pe)(i) 

3u 

P<1> 

■'4ep^ 

3? " 

REF 


uv 


(79) 


( 80 ) 


I- 


C-momentum 


,1/2 3? 


+ [(1 + 


1 i(i) 1 Iv ^1/2., 1 Cy3£ 


1/2 R 


1/2 


1/2 R'dK 


2 r (1) 3 V 1 1 3v 3 , o., 

^ H o 35 ■ 35^“’^ 

Rsref 


1 r r M 1 ^ 

+ ^{[v(l+^)l 35,5 

^®REF 


2 3u 7^ 91 I ,,(i) 

1/2 ^3^^ 3? ^ SraP ^ ^t’ 3c 


3?35 ^3T' ^3?^ 3C 


-2 ^ (1) usin0 ^ ^ 3n 

+ p<^> ^ + o(|H) (-^) |f> 

3C 

+ Ipu/^', - 1 „ (81) 


"®EEF 


+ (1 + — |) ^ ■35' ‘^^"^^“rEF *' 

^^REF 

((, +_J la . la , 

^ ^ l/2p^ 3? r^1/2i3 

R^rEF^ ^^REF^ 

C ^ y ^^REF,^.,^^ 3^ 

+ [ (1 + 1 / 2 ) Pr u Pr 2 

Pr 

. _Lrl_. CT , 1 Is .U (-, . £E N 1 ^ 
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Terms with superscript (i) are linearized terms and are updated after 
each iteration. 

Following reference 4, to enhance stability, the continuity 
equation is written for point B (figure 18) and the two momentum 
equations and the energy equation are written for point A. The 
resulting system of algebraic linear finite difference equations 
which is derived in detail in Appendix F is : 

in nH-l,n-l in m+l,n in m+l,n-M in m+l,n-l in nrl-l,n 


+ Fv,, +GT ,+HT +IT 

in nH-l,n+l in m+l,n-l in m+l,n in nrfl,n+l 


+ J, P , + K. p + L. p ,, 

in mi-l,n-l in m+l,n in m+l,n+l in 


(83) 


where : 


i = 1, 2, 3, 4 
2 < n < N-1 

A centered differencing has been used for the variables u, v, and T 

2 

in which the truncation error is of order [(A?) ], while for the 

pressure a forward two-point differencing has been applied. Though 

it is accurate only up to (A?) it gives a more stable solution than 

2 

a centered differencing which is accurate up to (A?) . This numerical 
instability associated with including the normal pressure gradient 
has been noted also by other investigators (reference 19) . 
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In matrix notation the system of equations (83) can be written 
in a more compact form: 

MX = i (84) 

where M is a block tridiagonal matrix defined: 


M = 



0 



0 0 
0 


\-2 


■2 \-2 


0 D M ^ 
N-1 N-1 


(85) 

Each element in M is a 4x4 matrix incorporating the coefficients of the 
system of the finite difference equations: 



®lk 


«ik 

“ik 


^Ik 

^Ik 

^Ik 

^Ik 

®2k 

®2k 

“2k 

^2k 


^2k 

^2k 

^2k 

^2k 

®3k 

^3k 

“3k 

Sk 

\ 

^3k 

“3k 

^3k 

^3k 

®4k 

^4k 

«4k 

^4k 


^4k 

^4k 

^4k 

^4k 
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and 


Ik 

’’ik 

“ik 

'^Ik 

2k 

“2k 

“2k 

*^2k 

3k 

“3k 

^3k 


4k 

“4k 

^^4k 

^^4k 


( 86 ) 


in equations (84) X is the unknown column vector. Each element of X 
is itself a 4-element column vector defined; 


u 

nrtljk 

^nrf-ljk 
T 

nrfljk 
^nrl-l,k 

(87) 

g in equation (84) is a known column vector whose elements are de- 
fined by: 



®k 


Ik 


2k 


3k 


4k 


( 88 ) 


Writing the equations for 2 < n < N-1 (see Appendix G) 4N-8 equations 

in 4N-1 unknowns are obtained, i.e., u ,, » v , T for 

m+l,n m+l,n m+l,n 
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1 < n < N and p for 2 < n < N. After applying the three boundary 

- — nrrl , n — — 

conditions at the wall (n=l) and the four conditions at the edge of 
the boundary layer (n=N) , an additional calculation for the pressure 
at the wall is needed in order to obtain the density at the wall 
through the equation of state. The induced pressure at the wall is 
calculated by applying the C~niomentum equation at the wall to obtain 
the pressure gradient and using a two-point differencing scheme the 
wall pressure can be determined and hence also the density at the 
wall. A detailed derivation of the expressions for calculating the 
pressure at the wall is given in Appendix H. 


Method of Solution 

Using the technique for block tridiagonal systems the block 
tridiagonal matrix M in equation (84) is decomposed into: 

M = LU (89) 


where L is the lower block diagonal matrix, and U is the upper block 
diagonal matrix. L and U are defined by: 


D 


L = 


0 0 


3 

0 D 


A^ 0 


4 \ 


^N-2 ^N-2 


0 

0 


DxT 1 

N-1 N-1 


0 


(90) 
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and 


U = 


2 2 


0 0 

N„ 


1 N 
■^N-3 N-3 


N, 


N-2 N-2 


^N-l 


(91) 


By comparison of the corresponding elements from both sides of 
equation (89) the following relationships are obtained: 

^2 ^ ^2 

^2 "" ^2 ^2 


(92) 


and; 


A = M - D N , 
n n n n-1 


3 < n < N-1 
2 < n < N-2 

The method for solving equation (84) is then to write it as 


N = A E 
n n n 


LUX = g 

or letting UX = W 
Equation (93) becomes: 
LW = i 


(93) 

(94) 


( 95 ) 
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The unknown vectors W are solved as follows: 



W = A ^ (i - D W J for 3 < n < N-1 (96) 

n n n n n-1 - - ' 

After solving for W^, the required unknown vectors can be obtained 

using the following procedure: 


Vi “ 


X 


n 


N-1 


= W -NX,, 
n n n+1 


N-2 > n > 2 (97) 

By underrelaxing the solution vectors the coefficients of the linear 
difference equations are updated and thus the iteration procedure 
is continued until convergence of the solution vector is obtained. 
The equation for the underrelaxation procedure is: 


X = X 
n n^ 


+ (X 


- X 


ru. 


‘OLD “NEW ’^OLD 

where is a positive number that is less than one. Most calculations 


) 


(98) 


were performed with 0,75. 


F. Numerical Results 

Numerical solutions using the second-order inner region pro- 
cedure have been obtained for a) the compressible laminar boundary 
layer on a flat plate over the Mach number range from 0 to 4, b) 
the compressible laminar boundary layer on a 20 degree half-angle 
cone at a free steam Mach number of 2.8, c) the compressible laminar 
boundary layer over the waisted body of reference 8 also at M^ = 2.8, 
and d) the turbulent boundary layer on a flat plate at M^ = 0. 

The laminar cases were computed for Pr = 0.72 and Sutherland viscosity. 
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To get the inner region solution started on its march down- 
stream, data are required along an initial station in addition to 
the specification of boundary conditions. In the present system of 
equations, profiles for u, v, T and p must be given. For the pre- 
sented laminar cases which all begin with nearly aero pressure 
gradient, the Blasius solution corrected for compressibility through 
the Howarth-Dorodintsyn transformation and adjusted when necessary 
for axial symmetry by the Mangier transformation has been used. The 
temperature profile has been obtained from the velocity profile 
using the Crocco integral and the induced pressure had been assumed 
to be zero across the boundary layer. Once the u, T and p profiles 
are known, the initial v profile is obtained by integration of the 
continuity equation. Sometimes, particularly at high Mach number, 
iteration of these approximate initial profiles is required in order 
to proceed downstream in a stable and convergent manner. 

a) Laminar Flat Plate ; The results obtained for the flat plate 
are in excellent agreement with classical boundary layer solutions. 
This is quite understandable since curvature effects are unimportant 
in this case. Nonetheless the present program is in physical co- 
ordinates and also gives direct calculation of the normal velocity 
distribution. In figures 19 and 20 the longitudinal and normal 
velocity distributions at M = 0 are seen to be essentially in- 
distinguishable from their Blasius counterparts. Similarly for the 
variation of the normal velocity at the edge of the boundary layer 
at M = 0 (figure 21) , except for a small blip at the beginning 
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indicating an initial datum line defect which however is damped out 
immediately. With increase in Mach number the calculation procedure 
yields convergent, smooth normal velocity distributions. The vari- 
ation of their edge values with distance along the plate at different 
Mach numbers is shown in figure 22. The induced pressure distribution 
is very small in this case as expected. 

The skin friction coefficient, form factor and displacement 
thickness distributions along the plate are shown in figure 23 for 
M = 0 for a range of unit Reynolds numbers and are seen to be 
indistinguishable from the similarity solutions. The very very 
slight differences are attributable to the different viscosity laws 
used, i.e. a linear viscosity temperature assumption in the Howarth- 
Dorodnitsyn stretching of the similarity solution as compared to the 
Sutherland viscosity relation in the finite difference solution. The 
same information at a Reynolds number of 1.5 x 10 but at 

Mach numbers up to 4 is shown in figure 24. These results again 
show good agreement with the similarity solutions. 

b) Cone : As a first axisymmetric example, numerical solutions 

for the compressible laminar boundary layer over a 20 degree half- 
angle cone, were obtained for a free stream Mach number of 2.8. The 
normal velocity at the edge of the boundary layer is given in figure 
25 which indicates that a stable and convergent solution has been 
achieved. As shown in the normal velocity has a negative sign for 
this cone angle indicating that the streamlines are directed into 
the surface whereas in the case of the flat plate or more slender cones 
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the edge streamlines are directed away from the surface. The skin 
friction coefficient along the surface is given in figure 26 in- 
dicating that a stable and convergent solution has been achieved. 

c) Laminar Flow over Waisted Body ; The compressible laminar 
boundary layer over the waisted body of revolution of Winter, Smith 
and Rotta (reference 8 ) has been calculated for M = 2.799. 

OQ 

Curvature effects, lateral as well as longitudinal are pronounced 
in this case. 

Variation of displacement thickness along the surface of = 2.799 
is shown in figure 27. The relatively large increase of the dis- 
placement thickness is mainly due to longitudinal curvature effect. 

In figure 28 the wall skin friction distribution is given and is seen 

to drop quite rapidly after X/h ^ 0.30. The normal velocity dis- 

REF 

tribution at the edge of the boundary layer is shown in figure 29. 

On the conical part of the waisted body the normal velocity is negative 
while further downstream the normal velocity at the edge changes 
sign due to curvature effects. 

The pressure distribution across the boundary layer is shown 

in figure 30 at three stations: “ 0.1125 on the conical part, 

X/L =0.325 near the maximum diameter and at X/L = 0.450 
REr REr 

after the inflection point of the surface. As shown the normal 
pressure gradient in the conical part is quite small whereas near 
the peak where curvature is pronounced the difference in wall press- 
ure can be of the order of 5% of the inviscid wall pressure. Down- 
stream at = 0.450 the pressure distribution across the boundary 
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layer changes due to curvature and differs noticeably from the 
inviscid value. 


d. Incompressible Turbulent Flow Over Flat Plate 

Numerical solutions have been obtained for the incompressible 
turbulent boundary layer over a flat plate using the two layer model 
for eddy viscosity. Initial profiles were generated using the law of 
the wall as formulated by Walz (reference 20) in terms of the follow- 
ing three algebraic relations corresponding respectively to the 
laminar sublayer for y+ < 4, a transition region for 4 < y+ < 26, 
and the logarithmic law for y+ > 26. 

u+ = y+ ^ 

— ay+ 

u+ = c^ • ln(l + y+) + c^ + [ (1 - “ C 2 *a)yH — ^^2^ ^ 

4 < y+ < 26 (100) 

u+ = c^£n y+ + y*” ^ (101) 

In these expressions 



y+ = 


u = A /p 


w 


The constants c^, c^ and a in equation (100) are taken as 2.50, 5.10 


and 0.3 respectively. 
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All the following results are for a reference Reynolds 
5 

number of 1.588 x 10 . In figure 31 are shown distributions of 

skin friction coefficient^ displacement and momentum thicknesses as 
well as edge normal velocity. The skin friction coefficience compare 

well with the Karman-Schoenherr relation as given in reference 21. 

Some oscillations are seen in the beginning of the edge normal velocity 
distribution but these rapidly die out. Obviously the initial 
datum line has some inconsistencies with the difference equations. 

Any such deviations in the initial data show up immediately in the 
normal velocity and induced pressure profiles. The displacement and 
momentum thickness distributions generally behave as expected. 

Figure 32 displays the u and v velocity profiles along the plate. 

The u-profiles that were calculated are as expected while for the 
normal velocity profiles a slight oscillatory character can be 
seen in the first three stations. 

G. Convergence and Stability 

The objective of this section is to list some of the major 
difficulties that have been met while trying to solve the non-linear 
partial differential equations that were replaced for solution by a 
system of linear algebraic equations using an implicit finite 
difference scheme. 

It seems that in a complex system of non-linear partial diff- 
erential equations such as developed herein for the inner layer, it 
is quite difficult to estimate stability criteria in closed form 
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by applying the VonNeumann method or some equivalent numerical sta- 
bility analysis. Also when the solution becomes unstable it is 
difficult to trace the cause for that instability. Therefore only 
by trying different grid sizing and different weighting coefficients 

A , X , X , and A , have regions of stable calculation been found, 
c u V T p 

The accuracy of the numerical solution has been checked with cases 
for which exact solutions exist. 

Continuity Equation: It has been recognized that the continuity 

equation written as a central differencing in the ? direction may 

lead to strong oscillations or even to divergence because of the 

boundary condition at the wall for the density or pressure (reference 

22) . One way to overcome this oscillatory behavior is to add an 

artificial eddy viscosity term into the continuity equation. In 

the present analysis however, a weighting factor of ” 0.85 gives 

well behaved results (see Appendix F) . Figures 33, 34, and 35 shows 

the normal velocity at the edge of the incompressible laminar 

boundary layer along a flat plate for A = 0.50; A^ = 0.85 and 

A =1.20 respectively. For A^ = 0.5, the solution is clearly 

oscillatory and with increasing amplitude of oscillation. The 

solution for A =1.20 (figure 35) is overdamped. A value of 
c 

A^ = 0.85 is nearly optimal since only slight overdamping is seen 
for the first two stations and further dcvmstream the results are 
the same as for the exact solution. The oscillatory behavior for 
A^ = 0.5 is also maintained for the compressible boundary layers. 

A representative result at = 0.5 is seen in figure 36. All inner 
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layer calculations other than these stability and convergence tests 

were performed with = 0.85. 

The solutions of Reyhner and Flugge-Lotz (reference 4) which 

utilize central differencing (X = 0.5) display the same oscillatory 

c 

behavior as seen in figures 33 and 36. For the purposes of inter- 
active calculation their edge normal velocity distribution was taken 
as the average of the peak-to-peak oscillations. They did not succeed 
in obtaining a smooth longitudinal variation of edge normal velocity. 


Pressure 

It has been recognized and reported (references 19 and 23) that 
departures from the proper solution can occur as a result of the 
(3p/3?) terms in the ? momentum and the energy equations. Besides 
in the present system there is an additional pressure gradient term (9p/3^) 
in the momentum equations and in the energy equation. It is found 
that these departures can be controlled by splitting the pressure 
as follows; 


p(?,?) = Pj.(5,?) 


( 102 ) 


where p (5) is the external pressure imposed on the viscous layer 
EXT 

by some outer region solution and p^(5,(;) is the induced pressure 
which is generally small compared to p„„„(C). 

Basically two options are used for calculating the pressure 
gradient in the 5 direction. They are: 

a) a centered difference for point A (see figure 18) 

- P^ )/A£ 


■^a,o = + (Pi 


m+l,n 


ra,n 


"nrfl 


3 ? 
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where: r- is a presumably known function and 

3^ I,m+l,n 

is the unknown induced pressure. 

b) The induced pressure gradient is taken as the gradient of 
average induced pressure across the boundary layer as follows: 


^ = 9_rl ( 


Pj (5,C)dC] 

o nrfl , n 

In this option the induced pressure gradient is a known function 
which depends on the updated induced pressure distribution across 
the boundary layer. Both of these options have been used successfully 
in the present calculations. There is no clear preference among 
them at this time. 

9Pi 

For each of the above-mentioned options for 
derivative of induced pressure was treated in several ways: 
a) a centered-dif f erence scheme accurate to order of 


2 

(A? ) and for which a stable solution could not be achieved. In 
n 

marching downstream a growing oscillatory behavior developed that 
could not be controlled. 

b) a backward difference scheme which includes two-point differ- 
encing accurate to order of this option stable 

solutions have been obtained. 

c) a forward difference-scheme which includes two-point differ- 
encing accurate to order of (A^^). Also here stable solutions 
have been obtained. The forward differencing is preferred since it 
is most convenient for matching when the inner region solution pro- 
cedure is used in interactive calculations. 


Grid Size : From experience with the program, it has been 


realized that careful selection of , and AC must be made in 

nt+1 n 

order to obtain stable and convergent solutions. The choice has 
upper bound as well as a lower bound. Roughly speaking taking 
A^^j^ of the order of 4 to 5 times ACj^ will give stable solutions 


CHAPTER IV 


INTERACTION PROCEDURE 

In the interactive mode, following the suggestion of Eerri and 
Dash (reference 6) , the inner and outer regions are matched at 
grid points where the Mach number is approximately 1.2 (figure 37). 
The matching procedure is however differently structured. 

In order that the computer programs for each region interact 
properly, the two programs have been restructured and written in an 
overlay mode. Namely the computer program is divided into segments, 
a main segment and other segments of lesser heirarchy connected to 
the main segment like branches of a tree. Thus when the main segment 
calls for a program in a particular segment of lesser hierarchy all 
other segments of equivalent lesser hierarchy are ignored and only 
this segment is in the operating mode for use. Overlay operation 
is needed mainly for three reasons: a) to reduce the need for 

storage capacity b) to eliminate the rewriting of the program 
because of the similar name variables in the two programs, and c) 
for convenience, while calculating one of the regions, the computer 
program for the other region is not needed, and therefore can be 
stored in a scratch file. 

A general flow diagram for the interactive mode is given in 
figure 38. Flow diagrams for the component outer region and inner 
region programs are shown in figures 39 and 40 respectively. 
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A. Initial Datum Line 

In order to start the numerical calculation in the interactive 
mode, an initial datum line is required which must be consistent in 
an interactive sense between the inner and outer regions. For con- 
venience the initial datum line has been fixed in a region where 
the cone flow solution is still valid. Accordingly the procedure 
for creating the initial datum line for axisymmetric flow is as 
follows; 

1. Calculate the “inviscid" flow field using the outer 
region program. 

2. Solve the inner region flow field using edge 
boundary conditions taken from the outer region 
solution of step 1. 

3. At a station sufficiently downstream of the apex such 
that the inner region solution is well behaved, 

the Mach number 1.2 location is determined within 

the inner region and the slope of the local 

streamline (tan0 = — ) is obtained. 

u 

4. At the station chosen in step 3, the outer region 
is recalculated using the Taylor-Maccoll equations 
with the slope of the streamline at the inter- 
face of the two regions (from step 3) taken as a 
boundary condition in place of the usual surface 
slope. Thus the shock angle and the flow field 
variables are changed slightly doe to the inter- 
action. 

5. The portion of the outer region profile from 
Mach number 1.2 to the edge of the boundary 
layer was scaled in order to get continuous pro- 
files of the flow variables for this portion 

of the outer region. This was done simply by 
multiplying the outer region profile by a 
dimensionless form of the inner profile. 

The procedure then gives complete information at the initial datum 


station. 
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B. Matching Procedure 

The matching between the inner and outer regions requires an 

iterative procedure. The matching condxtions are continuity of the 

flow variables u, v, T and p. Discontinuities in the derivatives with 

respect fo ? of the flow variables about the matching grid points are 

due to the use of different numerical procedures on either side of the 

2 

match points and to inconsistencies of order (6/L^gp) or higher in 

the systems of equations describing the two regions. The numerical 

results indicate that these discontinuities are very minor. Referring 

to figure 37, the matching grid points A and B are in the transonic 

regime where the Mach number is of the order of 1.2. 

Point B is chosen at the intersection of an extended streamline 

through point A and a normal to the wall at station (m+1) . The 

slope e„ of the streamline at point B (station m+1) is assumed 
B 

initially to be equal to the streamline slope 6^ at point A 
(station m) , where the flow variables are all known. The flow vari- 
ables u„, V,,, T„ and p„ can be determined from information at point 
B B B ij 

C (station m) using the boundary point procedure of the outer region 
program. 

The inner region can now be calculated for station mfl using 
the following converged data for point B as "edge condition": 

'"(iirf-l,N) 

"^(ntfl,N) 

P(m+1,N) 
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Using the normal velocity gradient at the edge of the inner region 
as determined through the continuity equation, a new normal velocity 
Vg^ for point B is obtained from the inner region solution. The 
new normal velocity is used to correct the slope at B to: 

-1 

0 = tan • — - (103) 

®NEW *^B 

with the new stream slope for point B the entire process is repeated 
until convergence is obtained for the flow variables at point B. 
Following this, the flow field in the outer region at station m+l 
can be solved using the outer region system thus completing the 
calculation for station m+1. 

C. Numerical Results 

The interactive program has been applied to two examples: a) 

a supersonic flow over a 20 degree half angle cone at = 2.80 and 

b) supersonic flow over the waisted body described in reference 8 also 

at M = 2.80. Both cases were for compressible laminar flow. 

00 


Cone : Numerical results for the wall static pressure are 

very slightly different than the inviscid flow field results as ex- 
pected for this weakly interactive ca^e. The difference of the shock 
angle between the interactive solution and the outer region solution 
is shown in figure 41. A slight oscillation is noticed at about 


b. Waisted Body : Variation of static pressure difference 

relative to the inviscid static pressure is shown in figure 42 for the 
interactive and inner region solutions . The inner region solution 
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shows the reduction in surface pressure due to longitudinal curvature 
effects. This reduction is restored by the displacement effects in 
the interaction solution. 

In Figure 43 the Mach number profile from the interactive 
calculation is compared with the separately calculated inviscid and 
boundary layer profiles. The interactive profile is reasonably 
continuous in slope at the matching line between the inner and 
outer regions . The interactive calculation indicates a larger skin 
friction than the boundary layer calculation and the approach to the 
inviscid Mach number distribution is slower than one would expect 
from comparison with the non- interactive boundary layer profile. 

It is suspected that this solution is not fully converged to the 
weak interaction solution. 

Nevertheless, the interactive program does work and the 
indication from these results is that more experience is required 
in starting procedure, choice of matching location and convergence 
criteria to make the program fully useful as a design tool. 


CHAPTER V 


SUMMARY 

A method has been developed for analysis of viscous- inviscid 
interactions in supersonic flows. The outer supersonic region of the 
flow field is represented by a two-family method of characteristics. 
In the present scheme, inclusion of viscous and tonductive terms in 
the formulation allows calculation of the supersonic portion of 
boundary layers and the tracing of wave reflection within those 
supersonic regions. It has been shown that shock wave patterns can 
be altered by inclusion of the correction terms. The inner region 
is handled by a second order boundary layer system that includes 
longitudinal and transverse curvature effects and a normal momentum 
equation. This nonlinear partial differential system of equations 
of parabolic type have been solved numerically in the physical plane 
by replacing them by a system of linear algebraic equations using a 
modified Crank-Nicolson implicit finite difference scheme. The 
nonlinearity of inner region equations was taken into account by an 
iterative procedure until the difference between the flow variables 
for successive iterations is as small as desired. Generally 4 to 10 
iterations are needed to get a convergent solution. This depends 
however on the accuracy of information on the initial datum line. 
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An interesting and useful byproduct of the present work is the 
experience gained in overcoming the difficulties of obtaining a stable 
and convergent solution for the coupled parabolic non-linear partial 
differential set in the inner region. Other investigators have 
encountered these difficulties even for less complicated systems of 
equations. They overcame the difficulties either by decoupling the 
equations, i.e. , solving the equations successively or by arbitrarily 
smoothing the oscillatory results. It has been recognized by other 
investigators that the terms involving the pressure gradient (3p/8C) 
lead to instabilities in the numerical results. Besides in the present 
procedure an additional pressure gradient (3p/3?) appearing in the 
C-momentum equation and energy equation is another source for causing 
numerical instabilities. It is found in the present work that the 
numerical instabilities can be controlled by splitting the pressure 
into induced pressure where 

"external" pressure is obtained from the outer layer solution pro- 
cedure. The ? derivative of the induced pressure was treated by a 
two-point backward or forward difference scheme that leads to a 
stable and convergent solution. A centered difference scheme would 
lead to an unstable solution. Also it has been recognized here that 
the continuity equation which is a first order differential 
equation when written as a central differencing in the Z, direction 
leads to strong oscillation and divergence of the numerical solution. 

In the present work a weighted differencing scheme has been used. 

The oscillatory behavior that was particularly pronounced in the normal 
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velocity has been damped out by choosing a proper weighting factor. 

Numerical results obtained herein for compressible laminar 
boundary layers agree very well with exact solutions. Numerical 
results for the waisted body where longitudinal curvature is important 
leads to pressure variation in the normal direction across the boundary 
layer. The wall static pressures differ by as much as 5% from f’pxx* 

In the interactive mode the matching between cue inner and the 
outer regions requires an iterative procedure. The normal velocity 
is the key parameter for matching the two solutions in the transonic 
region where the Mach number is of the order of 1.2. Thus it is 
important to get a stable and convergent inner solution for the normal 
velocity. It is interesting to note that the inner region solution 
converges at a faster rate in an interactive mode compared to the 
solution of the inner region in the non- interactive mode. The 
matching procedure works satisfactorily and numerical results for a 
20 degree half angle cone and for the waisted body of Winter, Smith 
and Rotta at M =2.8 have been obtained. 

OO 

Supersonic viscous- inviscid interactive analysis is an important 
area that has received only limited attention to date. The inter- 
active procedure developed- herein should be quite useful in dealing 
with supersonic flow fields where separation and strong-shock inter- 
action effects are absent. Further work is required in order to make 
interaction procedures such as the present one more encompassing and 


versatile. 
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APPENDIX A 

TRANSFORMATION TO STREAMLINE-NORMAL COORDINATES 

The non-dimensional system of equations (7) to (11) written in 
cartesian coordinates (x,y) is transformed to a curvilinear coordinate 
system (s,n) where s and n are respectively the distance along a 
streamline and the distance normal to the streamline. In this form 
it is easy and straightforward to identify and retain the transport 
effects such as viscous .shear and heat flux only normal to the stream- 
lines . 


Transformation Relation: The transformation relations from cartesian 

coordinates (x,y) to the curvilinear coordinates (s,n) are (see 
figure 1) : 


n . Q J 

— — = COS0 sin0 — 

9x 9S 3n 


3 . ^ 9 . ^9 

^ + =“30 -j;;; 


(A-1) 


Also the velocity components (u,v) in cartesian coordinates expressed 
in terms of the velocity vector magnitude q, and flow direction, 0, 


are: 


u = q COS0 
V = q sin0 


(A-2) 
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After substitution of relations (A-1) and (A-2), the governing 
equations (7) to (11) when simplified and expressed in curvilinear 
system (s,n) are; 

Continuity 

„ 8 / ^ cTn . „ 9 / On , . „ 3 / . „ o. 

COS0 T~(pq COS0 y ) - sxn0 7~Cpq cos0y ) + sin0 ”(pq sxnOy ) 

os 9n dp 

3 o 

+ COS0 sin0y ) = 0 (A-3) 

Expanding and rearranging the terms, the resulting equation is: 


9 / N , 90 apqsin0 

^(pq) + pq — = - 

Longitudinal momentum 

9 9 

pq cos0[cos0 — (qcos0) - sin0 — (qcos0) ] 
ds dn 


(A-4) 


9 9 

+ pq sin0[sin0 — (qcos0) + cos0 — (qcos0)] 
9s 9n 

= - — [cose - sins 




(A-5) 


9^9^ 2 

Considering an ordering procedure, g^'2'/g~2 is of the order of (6/L) , 

therefore terms of this order are neglected in the present analysis. 

The resulting expression is: 


^^1 Re 


1/1 

1 — COS0 ■r~[yy cos0 ~(qcos0)J 
REF y^ 


' 9 9 

+ sin0 ~[ysin0 — (qcos0)]} 

dn dn 

Expanding equation (A-5) and rearranging terms yields: 

pq |j(qcos0) = - -j— COS0 sin6 ^ 


(A-6) 


(A- 7) 




^^EF 
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Normal momentum 


9 3 

pq cos6[cos9 — (qsin6) - sin9 •r-(qsin9) ] 
'I o dn 


9 9 

+ pqsin9[sin9 ~(qsin9) + cos9 ~(qsin9)] 
a s dn 

[sin9 + cos 9 -^] + AQ 
d s 9n 2 




where : 
AQ 


2 - |jtpy®<=ose IjCqsine)] 

9 9 

+ sin9 — [ysin9 — (qsin9)]} 

3n 3n 


Expanding and rearranging yields 


pq |^(qsin9) = - ^ 




sin6 


it ■ to **^2 

^”ref 


(A-8) 


(A-9) 


(A- 10) 


Energy 

For the simplified thermodynamic assumptions used herein, namely 


Pr 


REF 


Pr 


- 1 - 


REF 


= 1 


y = k 


the energy equation becomes: 

9T 9T 9T 9X 

pqcos9(cos9 sin9 — ) + pqsinO (sin9 H cos9 ~) 

3s 3n 3s 9n 

- ^^ [qcos9 (cos9 - sin9 -^) + qsin6(sin9 + cos9 •^) ] 

Y dS 3n 3s 3n 


= AQ, 


(A-11) 


where : 
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— 1 rl „ 9 , cr _ 9T, . . 9 . 9T_, 


■-REF REF y" 

(Y-I)M^EF 9 9 2 

+ ■::: {y[cos0 — (qcos9) - sinO — (qsin9) ] 

9n 


Re 


■REF 

. 2„,9 


4- 2y [sin"^9 (-^(qcos9) )‘^ + cos^0 (-^(qsin0) ) ^^] } 

rl n ^-n 


^9n 


9n 


(A-12) 


Expansion and rearrangement yields; 


pq 


II ^ _( y-l ) l£ 


9s 


q -T^ + AQ, 
Y 9s ^3 


(A- 13) 


Upon multiplication of equation (A- 7) by dos0 and (A-10) by sin0 
and then adding the equations, the following is obtained: 


la + __i_ Ip. Q 

9s „2 9s ^1 


(A-14) 


ym: 


REF 


where : 

Ql = AQ^ COS0 + AQ^ sin9 (A-15) 

Multiplication of equation (A-7) by (-sin9) and (A-10) by cos9 and 
then adding the equations yields the following: 

2 36, 3£ _ . 


„2 


9n 


^EF 

where: 

= sin0 AQ^ - cos 9 AQ^ 
The energy equation becomes : 


pq 


II _ Xy- 1) „ l£ 


9s 


^ 9s 


(A-16) 


(A-17) 


(A- 18) 


where : 


Qj - AQ3/(Y-l)t4^ 
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Transport terms 

The correction terms Q^, and dud to viscous shear and heat flux 
which are considered as source terms in the method of characteristics 
solution of the system of equations (A-4) , (A-7) , (A-10) and (A- 13) 
can be simplified in the following way: 


term 


Q = COS0 AQ + sin6 AQ 

JL, d. ^ 

= COS0 [yy^cosO ~(qcos0)] 

9 3 

+ sin0 ■^[ysinB •^(qcos0)]} 


^^REF 


(A-19) 


Expansion and cancellation of terms reduces the expression for to: 


Q = — + OCOSQV l£ _ sin^0> 

^1 Re„„ 9n 9n y 9n ^^^9n y 9n 

REF 


(A- 20) 


Q2 - term 


= sin0’AQ^ - cos0*AQ2 
= I - "— {■— COS0 |---[yy^cos0 |-(qcos0)] 

+ sine -^[Msin0 —(qcose)} - -^[yy°cos0 |^Cqsin9)] 


'REF 


9 9 

+ sin0 — [ysinO — (qsin0)]} 
gn on 


(A- 21) 
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APPENDIX B 


METHOD OF CHARACTERISTICS 


The chs,racteristic equations are derived here for steady viscous 
compressible rotational two-dimensional or axisynunetric configurations 
in supersonic flow. The independent variables chosen are the static 
pressure, p and the streamline direction, 6. 

For convenience equations (12) to (16) are repeated here as 
(B-1) through (B-5). 


Continuity 

9 r \ , 98 npq . „ 

— (pq) + pq — = - sine 
as on y 

s-momentum 


, , 1 8p ^ 

pq ^ + ~T~ Is ‘ 


la. 

9s 


^\ef 


n- mo men turn 


2 90 . 1 9p 

pq 1 r — = -Q 

^ Ss ,.„2 3 n ^2 


Energy 


pq 


9T 

9s 


^^REF 

(y-i) 




Equation of state 
P = PT 


(B-1) 


(B-2) 


(B-3) 


(B-4) 


(b-5) 
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Equation (B-1) can be written as follows: 


p 9s q 3s 8n y 

Equation (B-2) upon dividing through by pq becomes 


ila = 

q 9 s 



1 1 
2 29 s 

'Y^REF 


The equation of state in differential form is: 


(B-6) 


(B-7) 


1 9p _ i (B-8) 

p9s p9s T9S 

and upon dividing through by pqT the energy equation becomes: 

2 

1 9T _ Q ^ (Y-1) Jt. (B-9) 

-p 3s P<lT 3 Y pi' 


Substitution of (B-7), (B-8) and (B-9) into (B-6) yields: 


p 3s pqT 


^ (y-1) 

Y 


pT 3s 


} + 



1 3^ 

2 2 
\ef 


1e 

3s 


^ _ osinB 

3n ~ y 


or 


1 9P 

2 2 3 s 


(M^-1) 


+ 


_90 

9n 


Osin0 


^^\ef 


pqT 


pq 


(B-10) 


In order to obtain the characteristics the following system con- 
sisting of equations (B-10), (B-3) and identities for the differentials 

dp and d6 must be solved: 


- — ^ 


(B-13) 
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Upon expanding equation (B-12) the equation of characteristics is 
obtained: 
dn 


ds 


= itanX 


(B-14) 


where: 


A = sin ^ ^ Mach angle. 

The compatibility condition is obtained by solving equation 
(B-13): 

i£ de 


+ p Psine ^ Q 

YP ~ sinAcosA y _2 pqT 3 


pq 


M__^sinA 

REF 


dx 


cosAcos (0±A) 


0 


L P J 15 ) 

The + and - signs correspond to C+ and C- characteristic lines res- 
pectively. 

The non-dimensional form of the entropy relation in terms of 
temperature, T and the pressure, p is: 

(y-1) T p 

where : 


S = 


R 


The viscous shear and heat flux lead to entropy changes along a 

streamline. Using the energy equation (B-4) : 

3T „ , (y- 1) 1 8p 

S'" T P 3s 


Thus : 

2 

dS .. Y 3T 1 8p _ Y , (y- 1) 1 3p-, 1 9p 

ds (Y“1) T 9 s p 9s (Y“1) pTq ^3 Y pT Ss"^ p 9s 


1 - 


or upon further reduction 

2 .,2 


dS ^ 


pq 3 


The variation of stagnation temperature, along a streamline 
as derived from the energy equation (B-4) and s-momentum (B-2) is: 




S + ^1 




APPENDIX C 


CONICAL FLOWFIELD 


The conical flowfield is derived for steady, isentropic, and 
irrotational flow with cylindrical symmetry about the x-axis 
(reference 24) . For the mathematical construction of the conical 
flow pattern, since there is no length scale, all flow variables 
depend only on the ray angle from the apex or on the ratio: 

(C-1) 

The differential equations used for this potential flow are: 

(C-2) 


X 


0 ) = 


^ (irrotational flow) 

9x 8y ^ 


and 


2 . 

V 3x 
c 




c c 


(C-3) 


(C-4) 


' '(Y+1) * 

where u is the velocity component in the x-direction, v is the 
velocity component in the y-direction, c is the local speed of sound 
and c^ is the critical sound speed. 

Expressing the equation in terms of m equations (C-2) and (C-3) 
become respectively: 

r + “ Is - ° 

d Oi oCO 


and 
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, 1 2, 3u ^ 9v ,2 2. Sv 2 

(c - u ) — - 2uv -(e - v)co-— +cv=0 

3o) 3w 3o) 


(C-6) 


Clearly this pair is equivalent to one equation of second order for 
one function only. Equation (C-6) assumes a particularly useful 
form when v is introduced as a function of u. Thus from (C-5): 


0 ) = - 


3v 

3u) 

9ti 

3w 


9v 

3u 


(C-7) 


Differentiating (C-7) with respect to to 
^ 3(0 Su ~ 3u^3u 3(0 


or 


3u 

3(0 


3% 


(C-8) 


Substituting (C-8) into (C-5) yields the relation: 


3 V 
3 (0 


3v 

3u 


(C-9) 


Introducing (C-7) , (C-8) and (C-9) into equation (C-6) leads to the 

following particularly simple form of the Taylor-Maccoll equation: 

2 — ^ 

^ Iv ^ ^ ^ _(u±v_|ul (C_io) 

3u 3u ^ 


Eliminating c“ by means of equation (C-4) yields; 

, 2 (^)(u+v|2)^ 

9u _ 2 . + (-^) _ 2 3u 


V 


3v 


'3u' 


.2 _ 




(c-11) 


Along the cone surface the flow has the direction of tne ray 
(jj = — traced by the cone in the x,y plane and thus: 

y 
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(G-12) 


The conditions to be satisfied along the conical shock (reference 24) 


are given by: 


Cl “ St) cos^3c„ + ~ 


Y+1 


SH q 


V = (q - u) cotg3„„ (C- 

O on. 

where q is the dimensionless free stream velocity and 3„„ is 
o bn 

the conical shock-angle. In addition the initial slope along the 
shock is given by: 


dv _ V 
3u u-q 


(C-15) 


APPENDIX D 


DIMENSIONAL ANALYSIS OF INNER REGION EQUATIONS 

The highly viscous flow in the inner region is assumed to be 
steady, compressible, laminar or turbulent, two-dimensional or 
axisyimnetric flow over an adiabatic or non-adiabatic surface. The 
gas is assumed to be perfect. The system of equations is that pre- 
sented by Maslen (reference 12) consisting of continuity, com- 
pressible Navier-Stokes and energy equations. The equations are put 
into curvilinear coordinates in which measured along the 

surface and C* is measured normal to the surface. The equations 
in dimensional form are: 


continuity 


^(p*u*r*‘^) + ~[p*v*r*'^(l +1^)] - 0 




g-momentum 


... ,. Bu* . , , C*. 9u* . p*u*v* , 3p* 

y r\ r • V I r\ Ti .JU rv JU 


1 


R* 9?* R* 9gj 

9 

R'V - u* 

2,9u* , 9g* 


a R*(R*+C*) 8^*f^*^* R*+;* 


+ 


r* 

2 9 


R* + V* 


9g* 




R*+5* 


V , 2u*a 

)] - — ) (u*sin0 +v*cos9 )sin6 

R* w w w 


I l'* 1$ + 2^(u*sln9^-H.*cosO^)] 

r* 


(D-2) 
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C*-momentum 
p*u* 


R* 8v* , * * IZ* _ p*u* i£* ^ 

R*+C* 3C* ^ ^ 3;* R*+^A 


3?* 






2y* 


13* a. * 


R*+C* ^ R*+?* 


-) + 


R* 


(R*+?*) ^*CJ 35* 


3 , . .a.3u* , 
iy*r* Cr — r + 




^35* 


R*+C* 


)] 


R* + V* 


cos6(u*sin0 + v*cos9 ) - — — ^[y*(— — )+ y* 

2 cos^<,u sina^ t- V coso^^; 2 r*+5* ^ ^ 3i;A 


(u*sin9 + v*cos0 ) 
w w •, 

+ ^ 


(D-3) 


energy 


Shi + „*„* ^ 92^ + V* ^-) 

(R*+?*) 35* ^ 3 5* XR*+5*) 35* 35* 


.a . 3h* 

T-X II* 

a ^ ^ 1 

^ -( ^ 


R* 


3 5* Pr* 

, 3h* 


R*r*''y* 

9 / . 

VT) — -A* / n v*i* \ / 


■ (R*+5*) 35* Pr*('R*+5*) 


+ 


2 R* + V* 

. + + 2 ( - ^ - - 

Pr*(R*+5*) ^ ^ 35* ^ R*+5* 


35* 


4- ^(u*sin9^ + v*cos6^)^ + + ' R^+^ > 

■ r* 

3u.^ 

- t( MttI + -S + ^(u*sin9__ + V*COS0 ))^] 


3u* 


R* ^ -u* 

Ml 


R*+5* 


35* 


w 


w 


(D-4) 


Equation, of state 
p* = p*RT* 


(D-5) 


In these equations, the angle 0 is the local slope of the axisymmetric 

w 

d-J-* 

body relative to the normal to the axis cos0 = ■5752 , and R* is the longi- 

w d5* 

tudinal rsdius of curvature taken as positive for convex surfaces (ref. 11 ). 
Dimensionless Equations 

Reexamining the equations of motion using dimensional analysis, 
a set of equations is obtained constituting a higher approximation 
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than the classical boundary layer for the inner region, since terms of 
order unity and of order (6/1^^^^) are retained. 

The non-dimensional variables are defined as follows: 


F, = 


P = 


- 






REF 


6 ’ ^ L, 


R* 


r* u* V* 

r = r ; u = ; v - 


REF 


REF 


u. 


REF 


V, 


REF 


'^REF^REF 


. n= P* - T . , _ 

I ’ P „ 5 T - , h - 


REF 


REF 


h* 

\ef 


y* 

y = — 

*^REF 


(D-6) 


The transformation relations are: 

8 ^ ^ _1 _8_ 

.8 ^ _8_ = 1 A 

8;* 8?8C* 68C (D-7) 


using the non-dimensional variables and the transformation relation, 

the dimensionless equations are then obtained. 

continuity 


L- 

9? 


[pur 


] + 


V. 


REF 


REF 


^REF ^ 
5 8? 


[pvr^ 


(1 + 


_6 

^REF 



(D-8) 


By setting 

V L 
REF REF ^ 

^EF*^ 

it is established that in conserving mass, the normal velocity is of 
the following mangitude compared to the longitudinal velocity: 




REF 

Thus equation (D-8) becomes: 


9C 


[pur''] + |-[pvr'' (1 + ^ ^ 


~] 


(D-9) 
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The remaining equations are: 
g-momentum 
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energy 
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a 6 2 

+ — (usine + vcos8 )) } 

" ^EF " 


(D-13) 


where : 
Re 


REF 


t> u L 
^ REF REF REF 

is the reference Reynolds number 
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and 


u c 
^REF p. 


Pr, 
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^F 


is the reference Prandtl nximber. 


In order that viscous forces be of same order as inertia forces; 
2 


.._1_ = (_A_) 

Pci '‘T '' 

REF REF 


(D-14) 


The non-dimensional equations of motion obtained by retaining 

terms of order unity and terms of order (6/L__„) are now presented. 

REr 

By using the Taylor series expansion for: 


1 + 


A 


= 1 - 


— * 1 + 


^REF ^ 
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equations (D-10) to (D-13) become: 
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g- momentum 
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energy 
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(D-19) 

The energy equation in terms of temperature using perfect gas relation 


dh = c dT 
P 


(D-20) 


becomes : 
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+ ^^"^^\ef Cl +7 -^ So "^L^RaS 


(D-21) 
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Rewriting equations (D-16) to (D-21) in terms of Re^^gj- using relation 
(D-14) yields: 


i 
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continuity 
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g- momentum 


pu 
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(D-24) 
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REF 


REF 


_ /•] + i S.N _1_ A_rlLE -^l-l I 

- t-L 1 /O D>' pr Orl-D^ -Yi-j ^ 


y 9T 


1/2 R^ a a?^Pr gc-i ^ . l/2p_ R 8S 
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(D-26) 
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APPENDIX E 


TURBUT.ENT INNER LAYER 


The dimensionless time- dependent second-order boundary layer 
equations are; 
continuity 


0 9p , 3 / . 9 r cf /-I I 1 C\ 1 

It + + 372 ^ 
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E~ momentum 
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t- moment urn 


1 r 9v , 3v . . 1 3v., 
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, V 2 r/, ■ 1 2 u 3u, 
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'REF ”“REF 

Derivation of the second— order turbulent boundary layer equations 


follows by using the Reynolds procedure of representing each quantity 
by the sum of its time average and its departure or fluctuation from 


the time average; namely 


u = u 4- u' 

V = V + V* 

T = T + T’ 


P = P + p' 

p = p + p' (E-5) 


where the barred quantities are time-averaged and the primed quantities 
are the fluctuations. Viscosity fluctuations are neglected herein 
as they do not contribute to the leading turbulent transport effects. 
The time average of any of the quantities in equation (E-5) is 


defined by: 



+ T 


f (t + r' )dr' 


(E-6) 
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Also by definition: 
r t + 
f'= - / ° 


f ’ (t + t’ ) dx’ =0 


(E- 7 ) 


where t is a time period large enough to give a stationary character- 
ization to the turbulence. 

The indicated time averaging is best carried out on transport 
forms of equations (E- 2 ) , (E- 3 ) and (E- 4 ) obtained by using the 
continuity equation (E- 1 ) . 
g-momentum 

Multiplying equation (E- 2 ) by r"^ and equation (E- 1 ) by u and 
adding the resulting equations yields: 
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4-momentum 


Multiplying equation (E- 1 ) by — 7-7:7 and equation (E- 3 ) by r^ 


Re 


and adding the equations yields: 
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energy 

Similarly following the same procedure multiply equation (E-4) 


by T and add them to obtain: 
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REF 


(E-10) 


The following substitutions are introduced into equations (E-1) , 
(E-8), (E-9) and (E-10): 
pu = pu + (pu) ' 

pv = pv + (pv) * (E-11) 

The time average of the above mentioned equations yields: 
continuity 
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^-momentum 
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An ordering procedure shows that the leading apparent transport terms 
are the correlations u'v’ and v'T’ . Within this ordering the equations 
reduce to: 
continuity 


3 r — a, . 3 r — ,, . 1 C. Ot 

■^[pur ] +^[pv(l +37^i)r 1 = 0 


(E-16) 


Re 


REF 


g-momentum 


4 . -~n + i -S-nIH + - 1 -- Puv + 9r 

pu + pv(l + R^34 1/2 R ^ 3? 

"^REF 


_ , 1 Jui_r j. __L_ In TTi 

^ g [yr (1 “ a 3?^^ „ 1/2 R^ ^ ^ 

^ ^^REF ^^REF 


1/2 R 3? 
^ ^^REF 


u 3 / 1 p u'v' 

FT7(yr ) 


Re 


1/2 R 
REF 


(E-17) 
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APPENDIX F 


COEFFICIENTS OF THE FINITE DIFFERENCE EQUATIONS 


A detailed derivation of the coefficients for the finite difference 

equations is herein given. The linearized equations (79) to (82) are 

written in difference form and have been multiplied by A? , , in 

nrri 

order that magnitude of the coefficients be less sensitive to step 
size. 

continuity 

The continuity equation is 

|^[r%u] + ‘-[pvr'^d + — ^ f] = 0 (F-1) 

^®REF 

Following the suggestion of Reyhner and Flugge-Lotz (reference 4) 
the continuity equation is written for point B (figure 18) as follows: 
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n 
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The C derivative has been taken as a weighted average with Xc as the 

weighting factor. Ac = -^ corresponds to a centered differencing 

scheme. Expanding equation (F-2) and multiplying by A5 yields: 
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Written in coefficient form: 


Au "IBu +Cu +Pv +Ev 
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the coefficients for the continuity equation (F-3) are: 
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g- momentum 

The linearized g-momentum is: 
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The difference equation for point A is written (figure 18) as follows; 
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Rearranging terms and multiplying all the equation by A?^^ yields: 
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Equation (F-8) in coefficient form becomes 
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(F-9) 
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(F-10) 


C-momentum 

The linearized ^-momentum equation is: 
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The difference equation written for point A (figure 18) is: 
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Re^„^ m+l Re„TXTx R^xx-rx-r, 

REF REF REF 

(i) , _ ,, 

_ 2 .3r X (i)x,., ^nrH,n+l ^nrH,n-l 
^ ^ V (1+K)AC^ 

(1-X ) (v - V ) . (i) 
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Expanding equation (F-12) and multiplying by ^ yields :. 
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Equation (F-13) in coefficient form is: 
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energy 

The linearized energy equation is: 
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The difference form of equation (F-16) for point A (figure 18) 


IS : 


(pu) 


(T .. - T ) „„ (i) >-, , _1__ e>,^"^ "''nrH.n"'M : 

(x) m4-l,n m,n j_ (1 \l2 r' " 


A? 


nrfl 




Re 


2 ; 


REF 


, . 2 (i)/^EXT^^^ , ^^mfln_ “ '• 

= (y-DMeef- (--d? ^ aT,” ^ 

7 1 C (i) ,^p^^m+l,n+l ^m+l,Z 

+ (y-l)4 ,(l+^f) 


^®REF 

+ (1-X ) i^sun^flsj’, 

p KAS 

n 


n 


Ill 


^ tj 1/2 R^ 
^®REF 


1 .. ^u^^iirfl.n+l ^mfl,n-l^ , ^u^^m,n4-l ’^m.n.-l\ 

_^1/2_J^ (1+K)A?^ (1+K)AC ■■ 

ReRgFR ti n 




(i) 


X„2 

[ ;7(t. 


K(1+K)AC 


o^-^1 4.1 + Kl^i 1 - (1+K)T ) 

2 nrl-l,n+l nrf-l,n-l m+l,n. 


n 


(1-A )2 

+ r(T + KT - (1+K)T .)] 

K(1+K)A(;^ m,n+l m,n-l m,n 

n 


Pr 

+ ■*■ ”T/2 +-^ P^))]> 

■ ^ ^^MF ^ 


(i) 


4.1 - 1..^1 1> 4.1-1 l) 

P T mrfl , n+1 nrH ,n-l , T m»n+l m,n-l , 

L /'14.V\A>- J 


(1+K)AC 


n 


(1+K)A? 


n (F-17) 


Expanding equation (F-17) and multiplying by Ag yields; 
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Equation (F-18) in coefficient form is: 
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APPENDIX G 

BOUNDARY CONDITIONS IN MATRIX FORM 


For solving the system of equations, the boundary conditions 
must be incorporated in difference form in order to have a closed 
system. As an example a detailed description is given here for the 
adiabatic case. Writing equation (84) for n = 2 the following is 
obtained: 


MjX2 + E3X3 = 83 
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forward differencing yields: 
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using three point 
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r = I 12 

12 

^12 " ®12 " °12^m+l,l 

are adjusted values of the matrix elements due to the three-point 

evaluation of wall temperature in the adiabatic case. Similarly for 

the other elements H' , I' , S’ (n = 2,3,4). 

n^ ^2 ^2 

for 2 < n < N-1 equation (84) can be written as: 


D X , + M X + E X^ , . = g (G-3) 

n n-1 n n n n+1 °n 

For n = N-1 equation (84) is: 

A-2 ^-A-1 "" ®N-1 

Incorporating the boundary conditions at the edge of the boundary layer 
n = N it can be shown that ^ and do not change and can 

be written as follows : 
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APPENDIX H 

CALCULATION OF THE PRESSURE AT THE WALL 


Some investigators that include the normal momentum equation 
tend to assume a zero pressure gradient in the C direction at the wall. 
Upon looking at the ^-momentum equation and satisfying the equation 
at the wall i.e., applying the boundary conditions it can be shown 
that the normal pressure gradient at the wall is not zero. Apply 
the boundary condition to equation (F-11) the following equation is 
obtained: 
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Solving for Re^{^ and condensing somewhat 
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Using two”polnt differencing the pressure at the wall is obtained by: 


p(nH-l,l) = p(ntfl,2) - -^| 
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WALL STATIC PRESSURE. P/PREF 
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Axial Distance, 


Figure 8(c) Wall Static Pressure Distribution for Waisted Body 
of Winter, Smith and Rotta (reference 8) 



Axial Distance, x/L 


REF 


Figure 8(d) Wall Static Pressure Distribution for Waisted 
Body of Winter, Smith and Rotta (reference 8) 
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Figure 9 


Streamline Patterns for the Waisted Body 


of Winter, 


Smith and Rotta (Ref. 8) at = 1.70 




Figure 11(a) Static Pressure Distribution Across the Flowfleld at - 1.70 
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Station Along Normal to Streamline 
<y-yB)/(y5H:yB) 
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Figure 14 Static Presaure Distribution Along the Cowl 


Surface at M 
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Figure 16 Coordinate System for Inner Region 
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Figure 19 Longitudinal Velocity Profile for Flat Plate 

Laminar Boundary Layer at = 0 
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Figure 20 Normal Velocity Profile for Flat Plate Laminar 
Boundary Layer at = 0 
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Figure 21 Normal ^/elo city Distribution at the Edge of Flat Plate 
Laminar Boundary Layer at M =0 
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Figure 23 Displacement Thickness, Form Factor and Skin 

Friction Distribution for Flat Plate Boundary 
Layer at = 0 
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Figure 26 Skin Friction Distribution of Compressible Laminar 
Boundary Layer Over 20 Degree Half Angle Cone at 
M,„ - 2.799 
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Figure 28 Skin Friction Distribution over Waisted Body 

at M = 2.799 
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Figure 29 Normal Velocity at the Edge of the Boundary 
Layer Over Waisted Body at = 


2.799 
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Figure 32(a) Longitudinal Velocity Profiles for Flat Plate Turbulent Boundary Layer at 

M 0 and = 1.588 x 10^ 
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Figure 32(b) Normal Velocity Profiles for Flat Plate Turbulent Boundary Layer at 

M •» 0 and Re„^_ = 1.588 x 10^ 
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Figure 33 Normal Velocity Distribution at the Edge of Flat Plate Laminar Boundary Layer 

at M sr 0 and X •*0.50 
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Figure 35 Normal Velocity Distribution at the Edge of Flat Plate Laminar Boundary 
Layer at M 0 and X ■■1.20 
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NORMAL EDGE l/ELOCITY, V/VREF 
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Figure 36 Normal Velocity Distribution at the Edge of Flat Plate Laminar Boundary 
Layer at = 0.50 and = 9.50 
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Figure 38 Flow Diagram for the Interactive Mode Computer Program 
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Figure 40 Flow Diagram for the Inner Region Computer Program 
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tAACH NOMBBR^ M 


Figure 43 Mach Number Distribution Across the Two Regions for the 

Waisted Body of Winter, Smith and Rotta (ref. 8) at 

M = 2.80 (y ^ . /y^ = 1.67): 
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